import pandas as pd  # 导入 pandas 库，用于处理 Excel 数据和表格操作
import numpy as np  # 导入 numpy 库，用于数值计算和数组操作
import gurobipy as gp  # 导入 gurobipy 库，用于构建和求解优化模型
from gurobipy import GRB  # 从 gurobipy 导入 GRB，用于访问优化状态和参数
import matplotlib.pyplot as plt  # 导入 matplotlib.pyplot，用于绘图
import matplotlib  # 导入 matplotlib，用于设置绘图参数
import time  # 导入 time 模块，用于记录程序运行时间

start_time = time.time()  # 记录程序开始运行的时间戳，用于后续计算运行时间
T = 24
# 设置随机种子以确保风电实际出力的可重复性
np.random.seed(0)  # 设置 numpy 的随机种子为 0，确保每次运行生成的随机数一致

# 加载 case39_data.xlsx 文件，包含 IEEE 39 节点系统数据
case39_file = "case39_data.xlsx"  # 定义 IEEE 39 节点系统数据的 Excel 文件路径
case39_data = pd.read_excel(case39_file, sheet_name=None)  # 读取 Excel 文件中的所有工作表，存储为字典

# 提取各个表格数据
bus_data = case39_data['bus']  # 从 case39_data 中提取 'bus' 工作表，包含节点数据
gen_data = case39_data['generator']  # 从 case39_data 中提取 'generator' 工作表，包含发电机数据
branch_data = case39_data['branch']  # 从 case39_data 中提取 'branch' 工作表，包含线路数据
gencost_data = case39_data['gencost']  # 从 case39_data 中提取 'gencost' 工作表，包含发电成本数据

# 提取 bus 表格的每一列为数组
bus_i = bus_data['bus_i'].values  # 提取 'bus_i' 列，存储节点编号
bus_type = bus_data['type'].values  # 提取 'type' 列，存储节点类型（1:负荷节点，2:发电节点，3:参考节点）
bus_Pd = bus_data['Pd'].values  # 提取 'Pd' 列，存储每个节点的负荷需求（MW）

# 提取 gen 表格的每一列为数组，包含所有 10 台发电机
gen_bus = gen_data['bus'].values  # 提取 'bus' 列，存储发电机节点编号
gen_mBase = gen_data['mBase'].values  # 提取 'mBase' 列，存储发电机基准容量（MVA）
gen_Pmax = gen_data['Pmax'].values  # 提取 'Pmax' 列，存储发电机最大出力（MW）
gen_Pmin = gen_data['Pmin'].values  # 提取 'Pmin' 列，存储发电机最小出力（MW）
gen_ramp_up = gen_data['ramp_up'].values  # 提取 'ramp_up' 列，存储发电机爬坡率上限（MW/h）
gen_ramp_down = gen_data['ramp_down'].values  # 提取 'ramp_down' 列，存储发电机爬坡率下限（MW/h）

# 提取 branch 表格的每一列为数组
branch_fbus = branch_data['fbus'].values  # 提取 'fbus' 列，存储线路起始节点编号
branch_tbus = branch_data['tbus'].values  # 提取 'tbus' 列，存储线路终止节点编号
branch_r = branch_data['r'].values  # 提取 'r' 列，存储线路电阻值
branch_x = branch_data['x'].values  # 提取 'x' 列，存储线路电抗值
branch_b = branch_data['b'].values  # 提取 'b' 列，存储线路电导值
branch_rateA = branch_data['rateA'].values  # 提取 'rateA' 列，存储线路功率容量限制（MW）

# 设置发电成本系数，与 DP 一致，包含所有 10 台发电机
gencost_a = gencost_data['cost_c2'].values  # 定义发电机的二次项成本系数（$/MW^2）
gencost_b = gencost_data['cost_c1'].values  # 定义发电机的一次项成本系数（$/MW）
gencost_c = gencost_data['cost_c0'].values  # 定义发电机的常数项成本（$）

# 读取系统负荷数据
load_file = "39节点24时刻系统负荷.xlsx"  # 定义系统负荷数据的 Excel 文件路径
df_load = pd.read_excel(load_file, header=0)  # 读取负荷数据，header=0 表示第一行为列名
Pdsum = df_load.sum(axis=1).values  # 计算每个时段的系统总负荷（MW），按行求和
Pd_nodes = df_load[[f'bus{i}_PL' for i in range(1, 40)]].values  # 提取每个节点的负荷数据（MW）

# 加载风电预测数据
wind_file1 = "39节点风电场1预测.xlsx"  # 定义风电场 1 预测数据的 Excel 文件路径
wind_file2 = "39节点风电场2预测.xlsx"  # 定义风电场 2 预测数据的 Excel 文件路径
wind_data1 = pd.read_excel(wind_file1, header=0)  # 读取风电场 1 预测数据
wind_data2 = pd.read_excel(wind_file2, header=0)  # 读取风电场 2 预测数据
wind_pre1 = wind_data1['pre_pow'].values  # 提取风电场 1 预测出力（MW），不强制切片
wind_pre2 = wind_data2['pre_pow'].values  # 提取风电场 2 预测出力（MW），不强制切片
# 验证风电数据长度
print(f"wind_pre1 数据长度: {len(wind_pre1)}")  # 调试：输出 wind_pre1 长度
print(f"wind_pre2 数据长度: {len(wind_pre2)}")  # 调试：输出 wind_pre2 长度
if len(wind_pre1) < 24 or len(wind_pre2) < 24:
    print(f"警告：风电预测数据不足 24 个时段！wind_pre1: {len(wind_pre1)}, wind_pre2: {len(wind_pre2)}")
    wind_pre1 = np.pad(wind_pre1, (0, 24 - len(wind_pre1)), mode='edge')  # 填充最后一个值
    wind_pre2 = np.pad(wind_pre2, (0, 24 - len(wind_pre2)), mode='edge')  # 填充最后一个值
wind_pre = wind_pre1 + wind_pre2  # 总风电预测出力
print(f"wind_pre 数据长度: {len(wind_pre)}")  # 调试：输出 wind_pre 长度



# 定义风电误差和概率，风电场 1 和 2 分别设置
C_reg_up = 120  # 定义上调备用成本（$/MW）
C_reg_down = 60  # 定义下调备用成本（$/MW）
wind_capacity1 = 198.5  # 定义风电场 1 最大容量（MW）
wind_capacity2 = 149.5  # 定义风电场 2 最大容量（MW）
wind_capacity = wind_capacity1 + wind_capacity2  # 总风电容量
wind_error_mean1 = 0  # 定义风电场 1 误差均值（MW）
wind_error_std1 = 10  # 定义风电场 1 误差标准差（MW）
wind_error_mean2 = 5  # 定义风电场 2 误差均值（MW）
wind_error_std2 = 10  # 定义风电场 2 误差标准差（MW）

# 计算灵敏度因子矩阵（Shift Factor Matrix）
NumBus = len(bus_i)  # 计算节点总数
NumBranch = len(branch_fbus)  # 计算线路总数
B_node = np.zeros((NumBus, NumBus), dtype=complex)  # 初始化节点导纳矩阵（复数矩阵）
B_line = np.zeros((NumBranch, NumBus), dtype=complex)  # 初始化线路-节点导纳矩阵（复数矩阵）

# 构建节点导纳矩阵和线路-节点导纳矩阵
for N in range(NumBranch):  # 遍历每条线路
    i = branch_fbus[N] - 1  # 获取线路起始节点索引（节点编号从 1 开始，转换为 0-based 索引）
    j = branch_tbus[N] - 1  # 获取线路终止节点索引
    Y = -1 / complex(0, branch_x[N])  # 计算线路导纳（仅考虑电抗，忽略电阻）
    B_node[i, i] += Y  # 更新节点导纳矩阵对角元素（起始节点）
    B_node[j, j] += Y  # 更新节点导纳矩阵对角元素（终止节点）
    B_node[i, j] += -Y  # 更新节点导纳矩阵非对角元素（起始-终止节点）
    B_node[j, i] = B_node[i, j]  # 确保节点导纳矩阵对称
    B_line[N, i] = -Y  # 更新线路-节点导纳矩阵（起始节点）
    B_line[N, j] = Y  # 更新线路-节点导纳矩阵（终止节点）
SF = np.real(B_line.dot(np.linalg.inv(B_node)))  # 计算灵敏度因子矩阵：线路-节点导纳矩阵乘以节点导纳矩阵的逆（取实部）

# 确定性求解火电机组出力
model = gp.Model("全时段确定性求解")  # 创建 Gurobi 优化模型，命名为“全时段确定性求解”
model.setParam('OutputFlag', 0)  # 减少求解日志输出

Ng = len(gen_bus)  # 定义发电机数量为 10
Pg = model.addVars(T, Ng, vtype=GRB.CONTINUOUS, name="Pg")  # 添加火电机组出力变量 Pg，维度为 T x n_gens，连续变量

# 设置火电机组出力上下限约束
for t in range(T):  # 遍历每个时段
    for i in range(Ng):  # 遍历每个发电机
        Pg[t, i].LB = gen_Pmin[i]  # 设置出力下限为发电机最小出力
        Pg[t, i].UB = gen_Pmax[i]  # 设置出力上限为发电机最大出力

# 定义目标函数：最小化总发电成本
cost = gp.quicksum(
    gencost_a[i] * Pg[t, i] * Pg[t, i] + gencost_b[i] * Pg[t, i] + gencost_c[i]
    for t in range(T) for i in range(Ng)
)  # 计算总成本：二次项 + 一次项 + 常数项，求和所有时段和发电机
model.setObjective(cost, GRB.MINIMIZE)  # 设置目标为最小化总成本

# 添加功率平衡约束
for t in range(T):  # 遍历每个时段
    model.addConstr(
        gp.quicksum(Pg[t, i] for i in range(Ng)) + wind_pre[t] == Pdsum[t])  # 约束：火电机组总出力 + 风电预测出力 = 系统总负荷

# 注释掉的爬坡率约束（未启用）
#for t in range(1, T):  # 遍历第 1 到 T-1 时段
#    for i in range(Ng):  # 遍历每个发电机
#        model.addConstr(Pg[t, i] - Pg[t - 1, i] <= gen_ramp_up[i])  # 约束：出力变化不超过上行爬坡率
#        model.addConstr(Pg[t - 1, i] - Pg[t, i] <= gen_ramp_down[i])  # 约束：出力变化不超过下行爬坡率

# 添加线路潮流约束
wind_bus1 = 6  # 风电场 1 连接节点 6
wind_bus2 = 16  # 风电场 2 连接节点 16
for t in range(T):  # 遍历每个时段
    for l in range(NumBranch):  # 遍历每条线路
        P_node = []  # 初始化节点注入功率列表
        gen_idx = 0  # 初始化发电机索引
        for n in range(NumBus):  # 遍历每个节点
            demand = Pd_nodes[t, n]  # 获取当前时段节点 n 的负荷
            node_power = -demand  # 节点注入功率初始化为负负荷
            for G in range(Ng):  # 遍历每个发电机
                if gen_bus[G] == n + 1:  # 如果发电机连接到当前节点
                    node_power += Pg[t, G]  # 加上发电机出力
                    gen_idx += 1  # 更新发电机索引
            if n + 1 == wind_bus1:  # 如果当前节点是风电场 1 节点
                node_power += wind_pre1[t]  # 加上风电场 1 预测出力
            if n + 1 == wind_bus2:  # 如果当前节点是风电场 2 节点
                node_power += wind_pre2[t]  # 加上风电场 2 预测出力
            P_node.append(node_power)  # 将节点注入功率添加到列表
        P_node = np.array(P_node)  # 转换为 numpy 数组
        flow = gp.quicksum(SF[l, n] * P_node[n] for n in range(NumBus))  # 计算线路 l 的潮流：灵敏度因子 × 节点注入功率
        model.addConstr(flow <= branch_rateA[l], name=f"Flow_Upper_t{t}_l{l}")  # 约束：线路潮流不超过上限（注释掉）
        model.addConstr(flow >= -branch_rateA[l], name=f"Flow_Lower_t{t}_l{l}")  # 约束：线路潮流不低于下限（注释掉）

# 求解优化模型
model.optimize()  # 调用 Gurobi 求解器求解模型

# 检查优化结果
if model.status == GRB.OPTIMAL:  # 如果优化成功
    Pg_values = np.array([[Pg[t, i].X for i in range(Ng)] for t in range(T)])  # 提取火电机组出力结果，存储为 T x n_gens 数组
    print("确定性求解优化成功！火电机组出力计划(MW)：")  # 输出成功提示
    for t in range(T):  # 遍历每个时段
        print(f"决策时刻 {t}: {Pg_values[t]}")  # 输出该时段的火电机组出力
if model.status == GRB.INFEASIBLE:  # 如果模型不可行
    print("模型不可行，开始写入 IIS")  # 输出不可行提示
    model.computeIIS()  # 计算不可行子系统（Irreducible Inconsistent Subsystem）
    model.write("model.ilp")  # 将不可行子系统写入文件 model.ilp
    print("已将不可行子系统写入 model.ilp 文件")  # 输出写入完成提示
    exit()

# 定义 MDP 的状态空间
state_space = []  # 初始化状态空间列表
for t in range(T):  # 遍历每个时段
    state = {  # 定义当前时段的状态
        't': t,  # 当前时段编号
        'Pd_t+1': Pdsum[t],  # 下一时段的总负荷
        'Pw_pre_t+1': wind_pre[t],  # 下一时段的风电预测出力
        'Pg_t': Pg_values[t]  # 当前时段的火电机组出力（从确定性求解结果获取）
    }
    state_space.append(state)  # 将状态添加到状态空间

# 定义动作空间
action_space = []  # 初始化动作空间列表
for t in range(T):  # 遍历每个时段
    action = {}  # 初始化当前时段的动作
    Pg_plan = np.zeros(Ng)  # 初始化火电机组计划出力数组
    for i in range(Ng):  # 遍历每个发电机
        Pg_plan[i] = max(gen_Pmin[i], min(gen_Pmax[i], Pg_values[t][i]))  # 确保出力在上下限范围内
        if t > 0:  # 如果不是第一个时段
            Pg_prev = state_space[t - 1]['Pg_t'][i]  # 获取上一时段的出力
            Pg_plan[i] = min(Pg_plan[i], Pg_prev + gen_ramp_up[i])  # 确保满足上行爬坡率约束
            Pg_plan[i] = max(Pg_plan[i], Pg_prev - gen_ramp_down[i])  # 确保满足下行爬坡率约束
    action['Pg_plan'] = Pg_plan  # 将计划出力存储到动作中
    action_space.append(action)  # 将动作添加到动作空间

# 参数设置
N_segments = 40  # 定义每个发电机的出力分段数
gamma = 0.95  # 定义折扣因子为 0.95，用于动态规划中的未来成本折现
C_reg_up = 120  # 重新定义上调备用成本（$/MW）
C_reg_down = 60  # 重新定义下调备用成本（$/MW）

# 初始化值函数参数
theta = np.zeros(Ng)  # 初始化权重参数 theta，维度为 K（每个发电机一个权重）
b_tk0 = np.zeros((T, Ng))  # 初始化值函数常数项，维度为 T x Ng
k = np.zeros((T, Ng, N_segments))  # 初始化值函数斜率，维度为 T x K x N_segments
P_seg = np.zeros((T, Ng, N_segments))  # 初始化出力分段变量，维度为 T x K x N_segments

# 出力分段设置
P_segments = []  # 初始化出力分段边界列表
rho_k = []
for i in range(Ng):
    segment_bounds = np.linspace(0, gen_Pmax[i], N_segments + 1)  # 从 0 到 gen_Pmax[k]
    segment_lengths = [(segment_bounds[n + 1] - segment_bounds[n]) for n in range(N_segments)]
    rho_k.append(segment_lengths)
rho_k = np.array(rho_k)
print(f"rho_k 分段长度: {rho_k}")
print("分段长度 rho_k：")
# 计算确定性求解的总成本
def compute_deterministic_cost():
    total_cost = 0  # 初始化总成本
    for t in range(T):  # 遍历每个时段
        for i in range(Ng):  # 遍历每个发电机
            Pg_k = Pg_values[t][i]  # 获取当前时段发电机 i 的出力
            cost_k = gencost_a[i] * Pg_k ** 2 + gencost_b[i] * Pg_k + gencost_c[i]  # 计算发电成本
            total_cost += cost_k  # 累加到总成本
    return total_cost  # 返回总成本

# 修复1的 calculate_theta 函数
def calculate_theta():
    global theta  # 声明 theta 为全局变量
    deterministic_costs = np.zeros((T, Ng))  # 初始化确定性成本矩阵，维度为 T x K
    for t in range(T):  # 遍历每个时段
        for i in range(Ng):  # 遍历每个发电机
            Pg_k = Pg_values[t][i]  # 获取当前时段发电机 i 的出力
            deterministic_costs[t, i] = gencost_a[i] * Pg_k ** 2 + gencost_b[i] * Pg_k + gencost_c[i]  # 计算发电成本
    total_cost = sum(sum(deterministic_costs))  # 计算所有发电机和时段的总成本
    theta = np.array([sum(deterministic_costs[:, i]) / total_cost for i in range(Ng)])  # 初始化 theta 为每个发电机的成本占比
    max_iterations_theta = 50  # 定义最大迭代次数
    epsilon_theta = 1e-4  # 定义收敛阈值
    alpha = 0.1  # 初始化学习率

    for iteration in range(max_iterations_theta):  # 迭代优化 theta
        weighted_costs = np.zeros(Ng)  # 初始化加权成本数组
        for i in range(Ng):  # 遍历每个发电机
            weighted_costs[i] = theta[i] * sum(deterministic_costs[:, i])  # 计算加权成本
        target_cost = sum(weighted_costs) / Ng  # 计算目标成本（平均加权成本）
        delta = np.zeros(Ng)  # 初始化偏差数组
        normalization = np.zeros(Ng)  # 初始化归一化因子数组
        for i in range(Ng):  # 遍历每个发电机
            delta[i] = weighted_costs[i] - target_cost  # 计算加权成本与目标成本的偏差
            normalization[i] = sum(deterministic_costs[:, i]) + 1e-6  # 计算归一化因子（避免除以 0）
        alpha = alpha / (1 + iteration * 0.1)  # 更新学习率（随迭代衰减）
        theta_new = theta - alpha * delta / normalization  # 更新 theta
        theta_new = np.maximum(theta_new, 0)  # 确保 theta 非负
        theta_new = theta_new / np.sum(theta_new)  # 归一化 theta 使总和为 1
        if np.max(np.abs(theta_new - theta)) < epsilon_theta:  # 检查是否收敛
            print(f"theta 在第 {iteration} 次迭代收敛")  # 输出收敛信息
            break
        theta = theta_new  # 更新 theta
    print(f"直接搜索法确定的权重 theta: {theta}")  # 输出最终 theta 值

# 初始化值函数
def initialize_value_function():
    global b_tk0, k, P_seg  # 声明全局变量
    V0_t = np.zeros(T)
    # 遍历每个时段 t，从 0 到 T-1
    for t in range(T):
        # 遍历未来时段 t_prime，从 t+1 到 T-1，计算从 t+1 到 T 的累积成本，文档公式24
        # 表示在决策后状态下，基于确定性调度结果的未来发电成本总和
        for t_prime in range(t + 1, T):  # 从 t+1 到 T
            for i in range(Ng):
                Pg_k = Pg_values[t_prime][i]
                cost_k = gencost_a[i] * Pg_k ** 2 + gencost_b[i] * Pg_k + gencost_c[i]
                V0_t[t] += cost_k
    for t in range(T):  # 遍历每个时段
        for i in range(Ng):  # 遍历每个发电机
            Pg_k = Pg_values[t][i]  # 获取当前时段发电机 i 的出力
            remaining_P = Pg_k  # 初始化剩余出力
            for n in range(N_segments):  # 遍历每个分段
                P_seg[t, i, n] = min(remaining_P, rho_k[i][n])  # 分配出力到当前分段
                remaining_P -= P_seg[t, i, n]  # 更新剩余出力
                if remaining_P <= 0:  # 如果剩余出力小于等于 0
                    break  # 退出分段分配

            base_slope = (2 * gencost_a[i] * Pg_k + gencost_b[i])  # 计算成本函数的斜率（导数）
            segment_idx = 0  # 初始化分段索引
            for n in range(N_segments):  # 遍历每个分段
                if sum(rho_k[i][:n + 1]) >= Pg_k:  # 找到出力所在的中心分段
                    segment_idx = n
                    break

            step_size = 0.1  # 定义斜率调整步长
            for n in range(N_segments):  # 遍历每个分段
                if n == segment_idx:  # 如果是中心分段
                    k[t, i, n] = -base_slope  # 设置斜率为负的基本斜率
                elif segment_idx > n >= 0:  # 如果是中心分段之前的分段
                    k[t, i, n] = -base_slope - (segment_idx - n) * step_size * base_slope  # 负斜率，递减
                else:  # 如果是中心分段之后的分段
                    k[t, i, n] = base_slope + (n - segment_idx) * step_size * base_slope  # 正斜率，递增

        if t == 0:  # 如果是第一个时段
            print(Pg_values[0][0])  # 输出发电机 1 在 t=0 的出力
            print(gencost_a[0], gencost_b[0])  # 输出发电机 1 的成本系数
            print(f"t=0 时刻火电机组 1 (i=0) 分段 1-20 的斜率 (k):")  # 输出提示
            slopes = [k[t, 0, n] for n in range(0, N_segments)]  # 提取发电机 1 在 t=0 的分段斜率
            print(f"  {slopes}")  # 输出斜率

    for t in range(T):  # 遍历每个时段
        for i in range(Ng):  # 遍历每个发电机
            slope_contribution = sum(k[t, i, n] * P_seg[t, i, n] for n in range(N_segments))  # 计算斜率贡献
            b_tk0[t, i] = V0_t[t] - slope_contribution  # 计算值函数常数项

    print(b_tk0)  # 输出值函数常数项矩阵

# 训练值函数
def train_value_function():
    global b_tk0, k, P_seg  # 声明全局变量
    initialize_value_function()  # 初始化值函数

    scenario_N = 1  # 定义训练场景数为 100
    lambda_val = 2  # 定义学习率参数 lambda
    beta = 1  # 定义学习率衰减参数 beta
    epsilon = 1e-4  # 定义收敛阈值
    alpha = lambda_val / (lambda_val + (scenario_N + 1) ** beta - 1)  # 计算初始学习率
    k_prev = k.copy()  # 复制当前斜率，用于后续收敛判断
    print(f"值函数训练参数：")
    print(f"  场景数目 (scenario_N): {scenario_N}")
    print(f"  步长参数 λ (lambda_val): {lambda_val}")
    print(f"  步长指数 β (beta): {beta}")
    for scenario_idx in range(scenario_N):  # 遍历每个场景
        scenario = []  # 初始化场景列表
        for t in range(T):  # 遍历每个时段
            wind_error1 = np.random.normal(wind_error_mean1, wind_error_std1)  # 风电场 1 误差
            wind_error2 = np.random.normal(wind_error_mean2, wind_error_std2)  # 风电场 2 误差
            Pw_actual1 = max(0, min(wind_capacity1, wind_pre1[t] + wind_error1))  # 风电场 1 实际出力
            Pw_actual2 = max(0, min(wind_capacity2, wind_pre2[t] + wind_error2))  # 风电场 2 实际出力
            Pw_actual = Pw_actual1 + Pw_actual2  # 总实际风电出力
            scenario.append(Pw_actual)  # 将实际风电出力添加到场景

        for t in range(T):  # 遍历每个时段
            Pg_t = state_space[t]['Pg_t']  # 获取当前时段的火电机组出力
            Pw_t = scenario[t]  # 获取当前场景的实际风电出力
            model = gp.Model("Value_Function_Update")  # 创建 Gurobi 模型，命名为“Value_Function_Update”
            model.setParam('OutputFlag', 0)  # 设置 Gurobi 不输出求解日志

            Pg_plan = model.addVars(Ng, vtype=GRB.CONTINUOUS, name="Pg_plan")  # 添加火电机组计划出力变量，维度为 K
            P_seg_vars = model.addVars(Ng, N_segments, vtype=GRB.CONTINUOUS, name="mu")  # 添加分段出力变量，维度为 K x N_segments
            reg_up = model.addVar(vtype=GRB.CONTINUOUS, name="reg_up", lb=0)  # 添加上调备用变量，非负
            reg_down = model.addVar(vtype=GRB.CONTINUOUS, name="reg_down", lb=0)  # 添加下调备用变量，非负

            for i in range(Ng):  # 遍历每个发电机
                Pg_plan[i].LB = gen_Pmin[i]  # 设置出力下限
                Pg_plan[i].UB = gen_Pmax[i]  # 设置出力上限

            if t > 0:  # 如果不是第一个时段（注释掉）
                Pg_prev = state_space[t - 1]['Pg_t']  # 获取上一时段的出力
                for i in range(Ng):  # 遍历每个发电机
                    model.addConstr(Pg_plan[i] - Pg_prev[i] <= gen_ramp_up[i])  # 约束：出力变化不超过上行爬坡率
                    model.addConstr(Pg_prev[i] - Pg_plan[i] <= gen_ramp_down[i])  # 约束：出力变化不超过下行爬坡率

            P_total = gp.quicksum(Pg_plan[i] for i in range(Ng))  # 计算火电机组总出力
            model.addConstr(
                P_total + Pw_t == state_space[t]['Pd_t+1']
            )  # 约束：火电机组总出力 + 实际风电出力 = 下一时段总负荷
            model.addConstr(reg_up >= state_space[t]['Pd_t+1'] - Pw_t - P_total)  # 约束：上调备用量
            model.addConstr(reg_down >= P_total - (state_space[t]['Pd_t+1'] - Pw_t))  # 约束：下调备用量

            for i in range(Ng):  # 遍历每个发电机
                model.addConstr(
                    gp.quicksum(P_seg_vars[i, n] for n in range(N_segments)) == Pg_plan[i]
                )  # 约束：分段出力总和等于计划出力
                for n in range(N_segments):  # 遍历每个分段
                    model.addConstr(P_seg_vars[i, n] >= 0)  # 约束：分段出力非负
                    model.addConstr(P_seg_vars[i, n] <= rho_k[i][n])  # 约束：分段出力不超过分段长度

            instant_cost = gp.quicksum(
                gencost_a[i] * Pg_plan[i] * Pg_plan[i] + gencost_b[i] * Pg_plan[i] + gencost_c[i]
                for i in range(Ng)
            ) + C_reg_up * reg_up + C_reg_down * reg_down  # 计算即时成本：发电成本 + 备用成本

            if t < T - 1:  # 如果不是最后一个时段
                value_function = gp.quicksum(
                    theta[i] * (
                            b_tk0[t + 1, i] + gp.quicksum(
                        k[t + 1, i, n] * P_seg_vars[i, n]
                        for n in range(N_segments)
                        )
                    )
                    for i in range(Ng)
                )  # 计算下一时段的值函数
                model.setObjective(instant_cost + gamma * value_function, GRB.MINIMIZE)  # 设置目标：最小化即时成本 + 折扣后的未来值函数
            else:  # 如果是最后一个时段
                model.setObjective(instant_cost, GRB.MINIMIZE)  # 设置目标：仅最小化即时成本

            model.optimize()  # 求解模型

            if model.status == GRB.OPTIMAL:  # 如果优化成功
                for i in range(Ng):  # 遍历每个发电机
                    for n in range(N_segments):  # 遍历每个分段
                        P_seg[t, i, n] = P_seg_vars[i, n].X  # 更新分段出力变量

                cost_observed = instant_cost.getValue()  # 获取即时成本
                if t < T - 1:  # 如果不是最后一个时段
                    next_value = sum(
                        theta[i] * (
                                b_tk0[t + 1, i] + sum(
                            k[t + 1, i, n] * P_seg[t + 1, i, n]
                            for n in range(N_segments)
                            )
                        )
                        for i in range(Ng)
                    )  # 计算下一时段的值函数
                    sample_value = cost_observed + gamma * next_value  # 计算样本值：即时成本 + 折扣后的未来值
                else:  # 如果是最后一个时段
                    sample_value = cost_observed  # 样本值仅为即时成本

                for i in range(Ng):  # 遍历每个发电机
                    for n in range(N_segments):  # 遍历每个分段
                        current_value = (b_tk0[t, i]
                                         + sum(
                                    k[t, i, m] * P_seg[t, i, m]
                                    for m in range(N_segments)
                                    if not np.isnan(k[t, i, m]) and not np.isnan(P_seg[t, i, m])
                                        ))  # 计算当前值函数
                        if P_seg[t, i, n] > 1e-3:  # 如果分段出力足够大（避免除以 0）
                            sample_slope = (sample_value - current_value) / P_seg[t, i, n]  # 计算样本斜率

                            k[t, i, n] = (1 - alpha) * k[t, i, n] + alpha * sample_slope  # 更新斜率（指数平滑）
                        if np.isnan(k[t, i, n]) or np.isinf(k[t, i, n]):  # 检查斜率是否无效
                            print(f"警告：无效 k 在 t={t}, i={i}, n={n}")  # 输出警告

        total_scenarios = scenario_idx + 1  # 更新已处理的场景数
        if total_scenarios % 1 == 0:  # 每处理一个场景
            print(f"扫描第 {total_scenarios} 个场景，斜率在t=0: {np.round(k[0, :, 0], 2)}")  # 输出 t=0 的斜率

        slope_diff = np.max(np.abs(k - k_prev))  # 计算斜率变化
        if slope_diff < epsilon:  # 如果斜率变化小于阈值
            print(f"在扫描 {total_scenarios} 个场景后收敛")  # 输出收敛信息
            break
        k_prev = k.copy()  # 更新上一轮斜率

    if total_scenarios == scenario_N:  # 如果达到最大场景数
        print(f"达到最大场景数 ({scenario_N})，训练停止")  # 输出停止信息

# 前向求解
def forward_scheduling():
    schedule = []  # 初始化调度计划列表
    total_cost = 0  # 初始化总成本
    for t in range(T):  # 遍历每个时段
        model = gp.Model("Forward_Scheduling")  # 创建 Gurobi 模型，命名为“Forward_Scheduling”
        model.setParam('OutputFlag', 0)  # 设置 Gurobi 不输出求解日志
        Pg_plan = model.addVars(Ng, vtype=GRB.CONTINUOUS, name="Pg_plan")  # 添加火电机组计划出力变量
        P_seg_vars = model.addVars(Ng, N_segments, vtype=GRB.CONTINUOUS, name="mu")  # 添加分段出力变量
        for i in range(Ng):  # 遍历每个发电机
            Pg_plan[i].LB = gen_Pmin[i]  # 设置出力下限
            Pg_plan[i].UB = gen_Pmax[i]  # 设置出力上限

        if t > 0:  # 如果不是第一个时段（注释掉）
            Pg_prev = schedule[t - 1]['Pg_plan']  # 获取上一时段的计划出力
            for i in range(Ng):  # 遍历每个发电机
                model.addConstr(Pg_plan[i] - Pg_prev[i] <= gen_ramp_up[i])  # 约束：出力变化不超过上行爬坡率
                model.addConstr(Pg_prev[i] - Pg_plan[i] <= gen_ramp_down[i])  # 约束：出力变化不超过下行爬坡率

        model.addConstr(
            gp.quicksum(Pg_plan[i] for i in range(Ng)) + state_space[t]['Pw_pre_t+1'] == state_space[t]['Pd_t+1']
        )  # 约束：火电机组总出力 + 风电预测出力 = 下一时段总负荷

        for l in range(NumBranch):  # 遍历每条线路
            P_node = []  # 初始化节点注入功率列表
            gen_idx = 0  # 初始化发电机索引
            for n in range(NumBus):  # 遍历每个节点
                node_power = -Pd_nodes[t, n]  # 初始化节点注入功率为负负荷
                for G in range(Ng):  # 遍历每个发电机
                    if gen_bus[G] == n + 1:  # 如果发电机连接到当前节点
                        node_power += Pg_plan[G]  # 加上发电机出力
                        gen_idx += 1  # 更新发电机索引
                if n + 1 == wind_bus1:  # 如果当前节点是风电场 1 节点
                    node_power += state_space[t]['Pw_pre_t+1'] * (wind_capacity1 / wind_capacity)  # 按比例分配风电场 1 出力
                if n + 1 == wind_bus2:  # 如果当前节点是风电场 2 节点
                    node_power += state_space[t]['Pw_pre_t+1'] * (wind_capacity2 / wind_capacity)  # 按比例分配风电场 2 出力
                P_node.append(node_power)  # 将节点注入功率添加到列表
            P_node = np.array(P_node)  # 转换为 numpy 数组
            line_powflo = gp.quicksum(SF[l, n] * P_node[n] for n in range(NumBus))  # 计算线路潮流
            model.addConstr(line_powflo <= branch_rateA[l], name=f"Flow_Upper_t{t}_l{l}")  # 约束：线路潮流不超过上限（注释掉）
            model.addConstr(line_powflo >= -branch_rateA[l], name=f"Flow_Lower_t{t}_l{l}")  # 约束：线路潮流不低于下限（注释掉）

        for i in range(Ng):  # 遍历每个发电机
            model.addConstr(
                gp.quicksum(P_seg_vars[i, n] for n in range(N_segments)) == Pg_plan[i]
            )  # 约束：分段出力总和等于计划出力
            for n in range(N_segments):  # 遍历每个分段
                model.addConstr(P_seg_vars[i, n] >= 0)  # 约束：分段出力非负
                model.addConstr(P_seg_vars[i, n] <= rho_k[i][n])  # 约束：分段出力不超过分段长度

        P_total = gp.quicksum(Pg_plan[i] for i in range(Ng))  # 计算火电机组总出力
        instant_cost = gp.quicksum(
            gencost_a[i] * Pg_plan[i] * Pg_plan[i] + gencost_b[i] * Pg_plan[i] + gencost_c[i]
            for i in range(Ng)
        )  # 计算发电成本

        if t < T - 1:  # 如果不是最后一个时段
            value_function = gp.quicksum(
                theta[i] * (
                        b_tk0[t + 1, i] + gp.quicksum(
                    k[t + 1, i, n] * P_seg_vars[i, n]
                    for n in range(N_segments)
                    )
                )
                for i in range(Ng)
            )  # 计算下一时段的值函数
            model.setObjective(instant_cost + gamma * value_function, GRB.MINIMIZE)  # 设置目标：最小化即时成本 + 折扣后的未来值函数
        else:  # 如果是最后一个时段
            model.setObjective(instant_cost, GRB.MINIMIZE)  # 设置目标：仅最小化即时成本

        model.optimize()  # 求解模型

        if model.status == GRB.OPTIMAL:  # 如果优化成功
            Pg_plan_values = np.array([Pg_plan[i].X for i in range(Ng)])  # 提取计划出力
            P_node = []  # 初始化节点注入功率列表
            gen_idx = 0  # 初始化发电机索引
            for n in range(NumBus):  # 遍历每个节点
                node_power = -Pd_nodes[t, n]  # 初始化节点注入功率为负负荷
                for G in range(Ng):  # 遍历每个发电机
                    if gen_bus[G] == n + 1:  # 如果发电机连接到当前节点
                        node_power += Pg_plan_values[G]  # 加上发电机出力
                        gen_idx += 1  # 更新发电机索引
                if n + 1 == wind_bus1:  # 如果当前节点是风电场 1 节点
                    node_power += state_space[t]['Pw_pre_t+1'] * (wind_capacity1 / wind_capacity)  # 按比例分配风电场 1 出力
                if n + 1 == wind_bus2:  # 如果当前节点是风电场 2 节点
                    node_power += state_space[t]['Pw_pre_t+1'] * (wind_capacity2 / wind_capacity)  # 按比例分配风电场 2 出力
                P_node.append(node_power)  # 将节点注入功率添加到列表
            P_node = np.array(P_node)  # 转换为 numpy 数组
            for l in range(NumBranch):  # 遍历每条线路
                line_flow = np.sum(SF[l, :] * P_node)  # 计算线路潮流
                if line_flow > branch_rateA[l] + 1e-6 or line_flow < -branch_rateA[l] - 1e-6:  # 检查潮流是否超出限制
                    print(f"警告：时段 {t} 线路 {l+1} 潮流不满足约束！")  # 输出警告
                    print(f"  实际潮流: {line_flow:.2f} MW, 上限: {branch_rateA[l]:.2f} MW, 下限: {-branch_rateA[l]:.2f} MW")  # 输出潮流信息
                    print(f"  起始节点: {branch_fbus[l]}, 终止节点: {branch_tbus[l]}")  # 输出线路信息
            schedule.append({'t': t, 'Pg_plan': Pg_plan_values, 'Pw_pre': state_space[t]['Pw_pre_t+1']})  # 将计划添加到调度
            total_cost += instant_cost.getValue()  # 累加即时成本
        else:  # 如果优化失败
            print(f"时段 {t} 优化失败，状态: {model.status}")
            model.computeIIS()
            model.write(f"model_t{t}.ilp")
            print(f"已将生成时段 {t} 的不可行子系统写入 model_t{t}.ilp 文件")
            return None, None
    return schedule, total_cost  # 返回调度计划和总成本

# 检验火电机组容量
def check_generator_capacity():
    total_gen_capacity = sum(gen_Pmax)  # 火电机组总最大容量
    net_load = Pdsum - wind_pre  # 净负荷 = 总负荷 - 风电预测出力
    print("\n火电机组容量检验：")
    print(f"火电机组总最大容量: {total_gen_capacity:.2f} MW")
    for t in range(T):
        if net_load[t] > total_gen_capacity:
            print(f"警告：时段 {t} 净负荷 ({net_load[t]:.2f} MW) 超过火电机组总容量 ({total_gen_capacity:.2f} MW)！")
        else:
            print(f"时段 {t}: 净负荷 = {net_load[t]:.2f} MW，火电机组容量充足")

# 执行训练和前向求解
deterministic_cost = compute_deterministic_cost()  # 计算确定性求解的总成本
print(f"确定性求解的总成本: {deterministic_cost:.2f}")  # 输出确定性总成本
calculate_theta()  # 计算权重 theta
train_value_function()  # 训练值函数
schedule, total_cost = forward_scheduling()  # 进行前向求解

# 输出结果
if schedule is not None:  # 如果调度计划存在
    print("最优调度计划：")  # 输出提示
    for item in schedule:  # 遍历每个时段的调度计划
        pg_vals = np.atleast_1d(item['Pg_plan'])  # 确保火电机组出力为一维数组
        pw_vals = np.atleast_1d(item['Pw_pre'])  # 确保风电预测出力为一维数组
        pg_str = ", ".join(f"{x:.2f}" for x in pg_vals)  # 格式化火电机组出力字符串
        pw_str = ", ".join(f"{x:.2f}" for x in pw_vals)  # 格式化风电预测出力字符串
        print(f"时段 {item['t']}: 火电机组出力（MW） = [{pg_str}]")  # 输出调度计划
    print(f"ADP求解的总成本（美元）: {total_cost:.2f}")  # 输出 ADP 总成本
else:  # 如果调度计划为空
    print("前向求解失败，无法生成调度计划")  # 输出失败信息

# 检验火电机组容量

# 绘制调度计划和净负荷曲线
matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 设置 matplotlib 使用 SimHei 字体，支持中文
matplotlib.rcParams['font.family'] = 'sans-serif'  # 设置字体族为 sans-serif
matplotlib.rcParams['axes.unicode_minus'] = False  # 确保负号正确显示

if schedule is not None:  # 如果调度计划存在
    time_steps = list(range(1, T + 1))  # 创建时间步列表 [1, 2, ..., T]
    Pg_plan_values = np.array([item['Pg_plan'] for item in schedule])  # 提取所有时段的火电机组计划出力
    Pw_pre_values = np.array([item['Pw_pre'] for item in schedule])  # 提取所有时段的风电预测出力


    plt.figure(figsize=(12, 8))  # 创建画布，设置大小为 12x8
    bar_width = 0.35  # 定义柱状图的宽度
    x = np.arange(len(time_steps))  # 创建时间步索引数组
    bottom = np.zeros(len(time_steps))  # 初始化柱状图底部高度
    for i in range(Ng):  # 遍历每个发电机
        plt.bar(x, Pg_plan_values[:, i], bar_width, bottom=bottom, label=f'火电机组 {i + 1}', alpha=0.7)  # 绘制火电机组出力的柱状图
        bottom += Pg_plan_values[:, i]  # 更新底部高度

    # 绘制风电预测出力曲线
    plt.plot(x, Pw_pre_values, 'g-', linewidth=2, label='风电预测出力', marker='o')
    # 绘制净负荷曲线
    plt.xlabel('T/h', fontsize=16)  # 设置 x 轴标签为“T/h”，字体大小 16
    plt.ylabel('功率/MW', fontsize=16)  # 设置 y 轴标签为“功率/MW”，字体大小 16
    plt.title('ADP算法调度结果', pad=20, fontsize=16)  # 设置标题，字体大小 16
    plt.xticks(x, time_steps)  # 设置 x 轴刻度为时间步
    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))  # 设置图例位置在画布外侧
    plt.tight_layout()  # 自动调整布局
    plt.savefig('adp_schedule_net_load.png')  # 保存图形
    plt.show()  # 显示图形

end_time = time.time()  # 记录程序结束时间
elapsed_time = end_time - start_time  # 计算程序运行时间
print(f"\n程序运行时间：{elapsed_time:.2f} 秒")  # 输出运行时间